Invariant /3-ensembles and the Gauss- Wigner crossover 
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We define a new diffusive matrix model converging towards the j3 -Dyson Brownian motion for all 
/3 £ [0, 2] that provides an explicit construction of /3-ensembles of random matrices that is invariant 
under the orthogonal/unitary group. For small values of /3, our process allows one to interpolate 
smoothly between the Gaussian distribution and the Wigner semi-circle. The interpolating limit 
distributions form a one parameter family that can be explicitly computed. This also allows us to 
compute the finite-size corrections to the semi-circle. 



Since Wigner's initial intuition that the statistical 
properties of the eigenvalues of random matrices should 
provide a good description of the excited states of com- 
plex nuclei, Random Matrix Theory has become one of 
the prominent field of research, at the boundary between 
atomic physics, solid state physics, statistical mechan- 
ics, statistics, probability theory and number theory [T]- 
3 . It is well known that the joint distribution of the 
eigenvalues of a large Gaussian random matrix can be 
expressed as the Boltzmann-Gibbs equilibrium weight of 
a one-dimensional repulsive Coulomb gas confined in an 
harmonic well. However, the effective "inverse tempera- 
ture" /? of the system cannot take arbitrary values but 
is quantized (in units of the repulsive Coulomb poten- 
tial). Depending on the symmetry of the random ma- 
trix, only three values are allowed (3 = 1 for symmetric 
real matrices, j3 = 2 for Hermitian matrices and (3 = 4 
for the symplectic ensemble. This is known as Dyson's 
"threefold way" . The existence of matrix ensembles that 
would lead to other, possibly continuous, values of [3, 
is a very natural question, and the quest for such en- 
sembles probably goes back to Dyson himself. Ten years 
ago, Dumitriu and Edelman j¥| have proposed an explicit 
construction of tri-diagonal matrices with non- identically 
distributed elements whose joint law of the eigenvalues is 
the one of /^-ensembles for general (3. Another construc- 
tion is proposed in p] p. 426-427] (see also [6]) and uses a 
bordering procedure to construct recursively a sequence 
of matrices with eigenvalues distributed as /3-ensembles. 
This construction gives not just the eigenvalue probabil- 
ity density of one matrix of the sequence but also the 
joint eigenvalue probability density of all matrices. This 
has lead to a renewed interest for those ensembles, that 
have connections with many problems, both in physics 
and in mathematics, see e.g. [21 [S]. The aim of the paper 
is to provide another construction of /3-ensembles that 
is, at least to our eyes, natural and transparent, and re- 
spects by construction the orthogonal/unitary symmetry 
[7j. Another motivation for our work comes from the re- 
cent development of free probability theory. "Freeness" 
for random matrices is the natural extension of indepen- 
dence for classical random variables. Very intuitively, 



two real symmetric matrices A, B are mutually free in 
the large N limit if the eigenbasis of B can be thought 
of as a random rotation of the eigenbasis of A (see e.g. 
[5] for an accessible introduction to freeness and for more 
rigorous statements) . "Free convolution" then allows one 
to compute the eigenvalue distribution of the sum A + B 
from the eigenvalue distribution of A and B, much in the 
same way as convolution allows one to compute the distri- 
bution of the sum of two independent random variables. 
In this context, the Wigner semi-circle distribution ap- 
pears as the limiting distribution for the sum of a large 
number of free random matrices, exactly as the Gaussian 
is the limiting distribution for the sum of a large number 
of iid (independent and identically distributed) random 
variables. A natural question, from this perspective, is 
whether one can build a natural framework that interpo- 
lates between these two limits. 

Let us first recall Dyson's Brownian motion construc- 
tion of the GOE [5] (for the sake of simplicity, we will 
only consider here extensions of the f3 — 1 ensemble, but 
similar considerations hold for (3 = 2 Hermitian matrices 
see [TU] for full details). One introduces a fictitious time 
t for the evolution of an N x N real symmetric matrix 
M(i). The evolution of the symmetric matrix is governed 
by the following stochastic differential equation (SDE): 

dM(t) = --M(t)dt + dH{t) (1) 

where dH(t) is a symmetric Brownian increment (i.e. 
a symmetric matrix whose entries above the diago- 
nal are independent Brownian increments with variance 
(dH.fj(t)) = ~(1 + 5ij)dt). Standard second order per- 
turbation theory allows one to write the evolution equa- 
tion for the eigenvalues A.; of the matrix M(i): 

1, , a 2 dt „ , . 

d\ l = --X l dt+ — V- r +<7d& i; (2) 

l l '—f. Ai — Aj 

where bi(t) are independent standard Brownian motions. 
This defines Dyson's Coulomb gas model, i.e. "charged" 
particles on a line, with positions Ai, interacting via a 
logarithmic potential, subject to some thermal noise and 
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confined by a harmonic potential. One can deduce from 
the above equation the Fokker-Planck equation for the 
joint density P({Xi},t), for which the stationary joint 
probability density function (pdf) is readily found to be: 
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with j3 = 1, independently of a 2 and where Z is a normal- 
ization factor. The above expression is the well known 
joint distribution of the eigenvalues of an TV x N ran- 
dom GOE matrix. The Wigner distribution can be re- 
covered either by a careful analysis of the mean marginal 
univariate distribution p(X) = J . . . J dA 2 . . . dA^P* (A = 
Ai, A 2 , . . . , Ajv) in the large N limit [12], or by using the 
above SDE ^ to derive a dynamical equation for the 
Sticltjes transform G(z,t) of p(A, t); 



G(z,t) 



1 N 

-Y 

N f-T 



-[ Ai(t) 



z e 



(4) 



With this scaling, the spectrum is spread out in a region 
of width of order y/~N and therefore z ~ a/A and G ~ 
1/y/N. Applying Ito's formula to G(z,t) and using ([2]), 
we obtain the following Burgers equation for G[T3]: 
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where a is introduced for later convenience, with a = 1 
for now. Note that we have neglected a noise term of 
mean zero and variance of order 1/NyN: G converges 
to the solution of Eq. ([5| only up to fluctuations of order 
(G 2 ) - (G) 2 ~ 1/A 3 (in agreement with e.g. Theorem 
9.2]). The contribution of these fluctuations is thus 1/N 
smaller than the diffusion term in Eq. ([5]), and can be 
neglected. 

For large N, the last (diffusion) term of Eq. ^ is of 
order 1/N smaller than the other ones. To leading order, 
the stationary solution (where the derivative with respect 
to time is set to 0) can be easily integrated with respect 
to z as: 
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aa 2 A^G 2 (z) + zG^z) = -1 . 



(6) 



where the integration constant comes from the boundary 
condition G(z) ~ —1/z when z — > oo. It is then easy 
to solve this equation to find the Stieltjes transform that 
indeed corresponds to the Wigner semi-circle density: 
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Now let us turn to the central idea of the present pa- 
per. In Dyson's construction, the extra Gaussian slice 



dM(t) that is added to H(f) is chosen to be indepen- 
dent of M(t) itself. The eigenbasis of dH(i) is a random 
rotation, taken uniformly over the orthogonal group. As 
mentioned above, this corresponds to free addition of ma- 
trices, and Eq. (|5| can indeed be derived (for N = oo) 
using free convolution [SJ. If instead we choose to add 
a random matrix dY(t) that is always diagonal in the 
same basis as that of M(t), the process becomes trivial. 
The diagonal elements of M(i) are all sums of iid ran- 
dom variables, and the eigenvalue distribution converges 
towards the Gaussian. The construction we propose is to 
alternate randomly the addition of a "free" slice and of a 
"commuting" slice. More precisely, our model is defined 
as follows: we divide time into small intervals of length 
1/n and for each interval [k/n; (k + l)/n], we choose inde- 
pendently Bernoulli random variables e£ , k £ N such that 
P[ £ £ = 1] = p = 1 - P[e£ = 0]. Then, setting e n t = ef nt] , 
our diffusive matrix process simply evolves as: 

dM n (t) = ~M n (t)dt + e?dH(t) + (1 - e?)dY(t) (8) 

where dH(t) is a symmetric Brownian increment as above 
and where dY(t) is a symmetric matrix that is co- 
diagonalizable with M„(t) (i.e. the two matrix have the 
same eigenvectors) but with a spectrum given by N in- 
dependent Brownian increments of variance cr 2 dt. It is 
clear that the eigenvalues of the matrix M„ (t) will cross 
at some points but only in intervals [k/n; (k + l)/n] for 
which e£ = (in the other intervals where they follow 
Dyson Brownian motion with parameter j3 = 1, it is well 
known that the repulsion is too strong and that colli- 
sions are avoided). In such a case, the eigenvalues are 
re-numbered at time t = (k + l)/n in increasing order. 

Now, using again standard perturbation theory, it is 
easy to derive the evolution of the eigenvalues of M n (t) 
denoted as AJ(t) < •• • < X%(t): 
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where the bi are independent Brownian motions also in- 
dependent of the eJJ, k £ N. 

A mathematically rigorous derivation provided in |10j 
allows one to show that the scaling limits Ai (t) , when n —> 
oo, of the eigenvalues A"(i) obey the following modified 
Dyson SDE: 



dt 
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with the additional ordering constraint Ai (t) < • • • < 
Ajy(i) for all t. One of the difficulty of the proof comes 
from the fact that when p < 1, there is a positive proba- 
bility for eigenvalues to collide in finite time (the ordering 
constraint is therefore useful at those points to re-start). 
The idea is then to show that collisions are in a sense 
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sufficiently rare for the above SDE to make sense (see 
[JDJ HI] for further details). Using the SDE (pi), one 



can derive as above the stationary distribution for the 
joint distribution of eigenvalues, which is still given by 
Eq. (|3| but with now /3 = a = p < 1. A very similar 
construction can be achieved in the GUE case, leading to 
/3 = 2p. As announced, our dynamical procedure, that al- 
ternates standard and free addition of random matrices, 
can lead to any /3-ensemble with f3 < 2. The correspond- 
ing matrices M(£) are furthermore invariant under the 
orthogonal (or unitary) group. This is intuitively clear, 
since both alternatives (adding a free slice or adding a 
commuting slice) respect this invariance, and lead to a 
Haar probability measure for the eigenvectors (i.e. uni- 
form over the orthogonal/unitary group). We have also 
proved that a collision leads to a complete randomization 
of the eigenvectors within the two-dimensional subspace 
corresponding to the colliding eigenvalues, see again [ID] . 

It is well known that the eigenvalue density corre- 
sponding to the measure P* given by ([!§ is the Wi gncr 
semi-circle for any /3 > 0. In fact, using (J5| with now 
a = (3 = p, one immediately finds that the eigenvalue 
density is a semi-circle with edges at ±as/2/3N. We sim- 
ulated numerically the matrix M n (i) with N = 200 for 
a very small step 1/n and until a large value of t so as 
to reach the stationary distribution for the eigenvalues. 
Then we started recording the spectrum and the nearest 
neighbor spacings (NNS) every 100 steps so as to sample 
the ensemble. We verified that the spectral density of 
M n (i = oo ) is indeed in very good agreement with the 
Wigner semi-circle distribution for f3 = 1/2. Our sample 
histogram for the NNS distribution is displayed in FigJI] 
We also added the corresponding Wigner surmise (which 
is expected to provide a good approximate description of 
the NNSD). 




FIG. 1. Empirical NNSD P(s) for the matrix M„(( = oo) for 
P = p — 1/2 with the Wigner surmise (red curve) correspond- 
ing to /3 = |, which behaves as when s — > 0. 



From the point of view of a cross-over between the 
standard Gaussian central limit theorem for random vari- 
ables and the Wigner central limit theorem for random 
matrices, we see that as soon as the probability p for a 
non-commuting slice is positive, the asymptotic density 



is the Wigner semi-circle, with a width of order \JpN . A 
continuous cross-over therefore takes place for p = 2c/ N 
with c strictly positive and independent of N. When 
c = 0, p(X) is a Gaussian of rms a, which indeed corre- 
sponds to the solution of Eq. ^ for a = 0. Setting for 
simplicity a = 1, the SDE for the system (Aj(i)) becomes 
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with the additional ordering constraint Xi(t) < ■ ■ • < 
Ajv(i) and the stationary joint pdf is still given by ([3]) 
but with now a vanishing repulsion coefficient /3 = 2c/N. 
In order to elicit the cross-over, we study Eq. Q with 
a = 2c/ N. The stationary differential equation corre- 
sponding to ^ (note this time that all terms are of the 
same and the second derivative term is not negligible) 
can be integrated with respect to z again as: 
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(12) 



where the integration constant comes from the boundary 
condition G ~ — 1/z for z — > oo. Note that (12 1 can be 



recovered directly from the saddle point equation route: 
under the measure P* with (3 = 2c/ N, the energy of a 
configuration of the Aj's can be expressed in term of the 
continuous state density p, neglecting terms <C 1, as: 
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A 2 p(A)dA-c / / ln(|A- A , |)p(A)p(A , )dAdA'. 



The probability density P* therefore rewrites in term of 
p as: 



P*[p] = Zexp —N 



£[p}+ pHp) 



S( p-1) 



where the entropy term, which is negligible when pi — p 
is of order 1, is now of the same order as the energy term 
(see [13] for a detailed discussion on the origin of the 
entropy term). We now need to minimize the quantity 
£ [p] + J pln(p) with respect to p. It is easy to see that 
the unique minimizer p c satisfies: 



Ap c (A) 



dA-2c 



PcWpc(x') 

(A-aO(A-A') 

P'cW 



d\d\' 



+ 



X - 



dX + v = 



where v is an integration constant. It is now straightfor- 



ward to derive ( 12 ) from this last equation by identifying 



each term and choosing the constant v so as to have the 
correct boundary condition for the Stieltjes transform of 
a probability measure. As expected physically, the dif- 



fusion term in (12 1 corresponds exactly to the entropy 



contribution to the saddle-point. 



Now, to solve (12), we first set G{z) := u'(z)/cu{z) to 
obtain a second order equation on u: 



i"(z) + zu'(z) + cu(z) = . 



(13) 
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It follows from the asymptotic behavior of G(z) that, for 
\z\ — >• oo, 



u(z) 
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(14) 



Eq. (131 can in turn be transformed with the change of 
function u(z) :— e~ z ^y(z) into a Schrodinger equation 
(sometimes called Weber differential equation) on y(z): 



y"(z) + [c 
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(15) 



The solutions of (15) are known (see [2lJ page 1067]) to 
write as y(z) — A 2 D c _i(z)+A 3 D_ c (iz) where D c _x,D_ c 
are parabolic cylinder functions and where A 2 and A 3 
are two constants. The general solution for u therefore is 



u(z) 



(A 2 D c -i(z) + A 3 D- C (iz)) and the correct 



asymptotic behavior of u is fulfilled for A 2 = 0. Now, 
one can recover the spectral density p c (A) associated to 
G by the classical inversion formula: 



pJX) = lim -S(G(A - ie)) 

1 i 



|S/(A)| 



(16) 
(17) 



where y% and y 2 are respectively the real and imaginary 



part of y. Using ( 15 ), it is straightforward to see that the 
derivative of (y 2 yi—y 2 y[) is identically zero and therefore 
the function is constant. This gives up to a normalizing 
constant the exact expression for p c (A). It happens that 
this constant can be computed explicitly |T5]. Eq. (12) 



was also studied in detail by Kerov [TB] and Askey & 
Wimp |T7j (see also [Hj). The final result for p(A) reads, 
for all c > 0: 



Pc(A) 
D- C (z) 



1 1 

\/2^r(i + c ) \D- C (i\)\ 2 

e -z 2 /4 /•"" 

r ' r Jo 



r(c) 



dxe zx " 2 x c 1 . 



(18) 
(19) 



This expression was again checked with numerical sim- 
ulations with very good agreement. The integral repre- 
sentation for D-c(z) does not hold for c = 0, but the 
function D- C {iu) is still well defined for all c € (— 1;0] 
(see [T7J Theorem 8.2.2]). It is easy to check that 
Po(u) — e~" / 2 /\/27r when c = 0, as expected. When 
c — > oo, the Wigner semi-circle law is recovered 
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At least when c is a non-negative integer, the integral 
form of D_ c (iA) can be computed analytically. This en- 
ables to find the tails of p c (u) for some values of c and 
large u. The asymptotic behavior reads: 
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(21) 




FIG. 2. Density p c (u) for c = 0,1,2,3,4, showing the pro- 
gressive deformation of the Gaussian towards Wigner's semi- 
circle. 



Let us return to fity for = a € (0;2). Interestingly, 
our method allows us to compute the correction to the 
Wigner semicircle for the spectral density for large but 
finite N due to the last diffusion term, which is usually 
neglected. Indeed one can solve as above the stationary 
equation of ^ keeping every term. This leads to the 
following corrected spectral density, valid for large but 
finite N: 
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(22) 



where a = 2/(2-/3) and c = f3N/(2 - 0). Numerical 
verification for this is provided in Fig. [3j 




FIG. 3. Numerical simulations of the state density of M n (t = 
oo) for p = /3 = 1/2, N — 50. Black curve: sample density. 
Red curve: Finite size correction ( 22 l , which coincides al- 
most perfectly with the sample density. Blue curve: N — > oo 
Wigner semi-circle density. 

The above discussion can also be formally extended 
to — 1 < c < 0, corresponding to a weakly attracting 
Coulomb gas (also mentioned in [7]; see also [TO] for an 
application). The system evolves, away from collisions, 
as in (11) but with c 6 (— 1;0] and with infinitesimal 
"up- down pushes" to separate particles when they col- 
lide. This kind of system has not been rigorously de- 
fined in the literature before but it should be possible, 
at least for a non trivial range of negative values of c, 
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since the attraction between particles is the same as for 
Bessel processes of dimension S G [0;1]. (see [20j page 
3]). We conjecture that the stationary density for large 
system is again given by the above Askey-Wimp-Kerov 
distributions p c but for the parameter range c G (—1; 0]. 
For c = — 1, the stationary density p_i is a Dirac mass 
at 0. Beyond this level, the attraction is too strong and 
the gas completely collapses on itself. 

As a conclusion, we have provided here the first explicit 
construction of invariant /3-ensembles of random matri- 
ces, for arbitrary (3 < 2. The stationary distribution for 
the eigenvectors is the Haar probability measure on the 
orthogonal group if < ft < 1, respectively unitary group 
if 1 < [3 < 2. We have found a natural scaling limit that 
allows one to interpolate smoothly between the Gaus- 
sian distribution, relevant for sums of independent ran- 
dom variables, and the Wigner semi-circle distribution, 
relevant for sums of free random matrices. The inter- 
polating limit distributions form a one parameter family 
that can be explicitly computed. Let us mention three 
interesting open problems. First, our alternate Bernoulli 
process of commuting and free matrix 'slices' can prob- 
ably be done differently, for example by introducing an 
Ornstcin-Uhlcnbcck process on the orthogonal group that 
mean-reverts towards the identity matrix. By using these 
matrices O to construct the 'slices' as dM(i) = O t dA O, 
one may be able to generate other interesting ensembles 
by tuning the parameters of the Ornstein-Uhlenbeck pro- 
cess, again interpolating between commuting and free ad- 
dition. Second, it would be interesting to know how the 
eigenvalue spacing distribution smoothly interpolates be- 
tween the Poisson distribution and Wigner's surmise. Fi- 
nally, the statistics of the largest eigenvalue is also very 
interesting (and now well known for j3 > 0, see [22Tl25j ): 
one should be able to interpolate smoothly, as a function 
of c, between the well-known Gumbel distribution of ex- 
treme value statistics and the Tracy- Widom(/3) distribu- 
tions. Whether this can be mapped into a generalized 
KPZ/Directed polymer problem remains to be seen. 
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